### Out of Sight, Out of Mind? ###
### Table 4 in the Main Text and its extended version as Table A15 in the Appendix ###

rm(list=ls())

## The package versions used to replicate the results are given next to each library.

library(magrittr) ## Version: 2.0.1
library(dplyr) ## Version: 1.0.5
library(readstata13) ## Version: 0.9.2
library(interflex) ## Version: 1.2.5
library(ggplot2) ## Version: 3.3.3
library(gridExtra) ## Version: 2.3

d <- read.dta13('./final_data.dta')


lograyic_models <- d %>% filter(year<2014) %>% select(rn, school_id_cluster, Dakp_3, Dwalk_3, lograyic09, Duniversity_3, Dodr_3, Dydr_3, Dpopulation_3, Dhospital_pri_3, Dhospital_pub_3)

lograyic_models <- lograyic_models[complete.cases(lograyic_models), ]


lograyic_models_ilgenel <- d %>% filter(year<2014) %>%select(rn, school_id_cluster, Dakp_5, Dwalk_5, lograyic09, Duniversity_5, Dodr_5, Dydr_5, Dpopulation_5, Dhospital_pri_5, Dhospital_pub_5)

lograyic_models_ilgenel <- lograyic_models_ilgenel[complete.cases(lograyic_models_ilgenel), ]

walk_lograyic_binning <- interflex(Y = "Dakp_3", D = "Dwalk_3", X = "lograyic09", Z = c('Duniversity_3', 'Dodr_3', 'Dydr_3', 'Dpopulation_3', 'Dhospital_pri_3', 'Dhospital_pub_3'),data = lograyic_models,estimator = 'binning', cl="rn", main = "", nbins = 3, Xlabel = "Property values (2009)", vcov.type = 'cluster', xlab = 'Moderator: Property values (2009)', ylab = 'Marginal effect of Dwalk on Dakp', cex.lab = 1, cex.axis = 1, theme.bw = TRUE, height=6, width = 4, ylim = c(-0.35, 0.2))


walk_lograyic_binning_ilgenel <- interflex(Y = "Dakp_5", D = "Dwalk_5", X = "lograyic09", Z = c('Duniversity_5', 'Dodr_5', 'Dydr_5', 'Dpopulation_5', 'Dhospital_pri_5', 'Dhospital_pub_5'),data = lograyic_models_ilgenel,estimator = 'binning', cl="rn", main = "", nbins = 3, Xlabel = "Property Prices (2009)", vcov.type = 'cluster')


## Gather both results under one table

# municipal council results
lograyic_binning_results <- summary(walk_lograyic_binning$model.binning)

lograyic_binning_results_estimates <- lograyic_binning_results$coefficients[ , 1]

lograyic_binning_results_se <- lograyic_binning_results$coefficients[ , 2]

# provincial council results
lograyic_binning_results_ilgenel <- summary(walk_lograyic_binning_ilgenel$model.binning)

lograyic_binning_results_ilgenel_estimates <- lograyic_binning_results_ilgenel$coefficients[ , 1]

lograyic_binning_results_ilgenel_se <- lograyic_binning_results_ilgenel$coefficients[ , 2]

## Table 4 ----
results_table <- data.frame(variables = c('Dwalk*low', '', 'Dwalk*medium', '', 'Dwalk*high', ''), m1 = c(lograyic_binning_results_estimates['D.G.1'] , lograyic_binning_results_se['D.G.1'] , lograyic_binning_results_estimates['D.G.2'],lograyic_binning_results_se['D.G.2'], lograyic_binning_results_estimates['D.G.3'],lograyic_binning_results_se['D.G.3']), m2 = c(lograyic_binning_results_ilgenel_estimates['D.G.1'] , lograyic_binning_results_ilgenel_se['D.G.1'] , lograyic_binning_results_ilgenel_estimates['D.G.2'],lograyic_binning_results_ilgenel_se['D.G.2'], lograyic_binning_results_ilgenel_estimates['D.G.3'],lograyic_binning_results_ilgenel_se['D.G.3']))

results_table$m1 <- round(results_table$m1, digits = 3)
results_table$m2 <- round(results_table$m2, digits = 3)

results_table[7, ] <- c('Observations', nrow(lograyic_models), nrow(lograyic_models_ilgenel))
results_table[8, ] <- c('Covariates', 'YES', 'YES')
results_table[9, ] <- c('Elec2009', 'MC', 'PC')

## Table 4
results_table

## Save Table 4 as .pdf
pdf("table4.pdf", height=11, width=8.5)
grid.table(results_table)
dev.off()


## Table A15: Expanded version of Table 4 -----

full_results_table <- data.frame(variables = c('PropertyLow', '', 'PropertyLow*PropertyValue2009', '', 'PropertyMedium', '', 'PropertyMedium*PropertyValue2009', '', 'PropertyHigh', '', 'PropertyHigh*PropertyValue2009','', 'Dwalk*PropertyLow', '', 'Dwalk*PropertyLow*PropertyValue2009', '', 'Dwalk*PropertyMedium', '', 'Dwalk*PropertyMedium*PropertyValue2009', '', 'Dwalk*PropertyHigh', '', 'Dwalk*PropertyHigh*PropertyValue2009','', 'Duniversity', '', 'DOldDependencyRatio', '', 'DYoungDependenyRatio', '', 'DPopulation', '', 'DHospitalPrivate', '', 'DHospitalPublic', ''), m1 = c(lograyic_binning_results_estimates['G.1'],lograyic_binning_results_se['G.1'],lograyic_binning_results_estimates['GX.1'], lograyic_binning_results_se['GX.1'], lograyic_binning_results_estimates['G.2'], lograyic_binning_results_se['G.2'], lograyic_binning_results_estimates['GX.2'],lograyic_binning_results_se['GX.2'],lograyic_binning_results_estimates['G.3'], lograyic_binning_results_se['G.3'], lograyic_binning_results_estimates['GX.3'], lograyic_binning_results_se['GX.3'], lograyic_binning_results_estimates['D.G.1'],lograyic_binning_results_se['D.G.1'], lograyic_binning_results_estimates['D.GX.1'], lograyic_binning_results_se['D.GX.1'], lograyic_binning_results_estimates['D.G.2'], lograyic_binning_results_se['D.G.2'], lograyic_binning_results_estimates['D.GX.2'],lograyic_binning_results_se['D.GX.2'], lograyic_binning_results_estimates['D.G.3'], lograyic_binning_results_se['D.G.3'], lograyic_binning_results_estimates['D.GX.3'], lograyic_binning_results_se['D.GX.3'],lograyic_binning_results_estimates['Duniversity_3'], lograyic_binning_results_se['Duniversity_3'], lograyic_binning_results_estimates['Dodr_3'], lograyic_binning_results_se['Dodr_3'], lograyic_binning_results_estimates['Dydr_3'], lograyic_binning_results_se['Dydr_3'],lograyic_binning_results_estimates['Dpopulation_3'],lograyic_binning_results_se['Dpopulation_3'], lograyic_binning_results_estimates['Dhospital_pri_3'], lograyic_binning_results_se['Dhospital_pri_3'],lograyic_binning_results_estimates['Dhospital_pub_3'], lograyic_binning_results_se['Dhospital_pub_3']), m2 = c(lograyic_binning_results_ilgenel_estimates['G.1'],lograyic_binning_results_ilgenel_se['G.1'],lograyic_binning_results_ilgenel_estimates['GX.1'], lograyic_binning_results_ilgenel_se['GX.1'], lograyic_binning_results_ilgenel_estimates['G.2'], lograyic_binning_results_ilgenel_se['G.2'], lograyic_binning_results_ilgenel_estimates['GX.2'],lograyic_binning_results_ilgenel_se['GX.2'],lograyic_binning_results_ilgenel_estimates['G.3'], lograyic_binning_results_ilgenel_se['G.3'], lograyic_binning_results_ilgenel_estimates['GX.3'], lograyic_binning_results_ilgenel_se['GX.3'], lograyic_binning_results_ilgenel_estimates['D.G.1'],lograyic_binning_results_ilgenel_se['D.G.1'], lograyic_binning_results_ilgenel_estimates['D.GX.1'], lograyic_binning_results_ilgenel_se['D.GX.1'], lograyic_binning_results_ilgenel_estimates['D.G.2'], lograyic_binning_results_ilgenel_se['D.G.2'], lograyic_binning_results_ilgenel_estimates['D.GX.2'],lograyic_binning_results_ilgenel_se['D.GX.2'], lograyic_binning_results_ilgenel_estimates['D.G.3'], lograyic_binning_results_ilgenel_se['D.G.3'], lograyic_binning_results_ilgenel_estimates['D.GX.3'], lograyic_binning_results_ilgenel_se['D.GX.3'],lograyic_binning_results_ilgenel_estimates['Duniversity_5'], lograyic_binning_results_ilgenel_se['Duniversity_5'], lograyic_binning_results_ilgenel_estimates['Dodr_5'], lograyic_binning_results_ilgenel_se['Dodr_5'], lograyic_binning_results_ilgenel_estimates['Dydr_5'], lograyic_binning_results_ilgenel_se['Dydr_5'],lograyic_binning_results_ilgenel_estimates['Dpopulation_5'],lograyic_binning_results_ilgenel_se['Dpopulation_5'], lograyic_binning_results_ilgenel_estimates['Dhospital_pri_5'], lograyic_binning_results_ilgenel_se['Dhospital_pri_5'],lograyic_binning_results_ilgenel_estimates['Dhospital_pub_5'], lograyic_binning_results_ilgenel_se['Dhospital_pub_5']))

full_results_table$m1 <- round(full_results_table$m1, digits = 3)
full_results_table$m2 <- round(full_results_table$m2, digits = 3)

full_results_table[37, ] <- c('Observations', nrow(lograyic_models), nrow(lograyic_models_ilgenel))
full_results_table[38, ] <- c('Elec2009', 'MC', 'PC')

## Table A15
full_results_table

## Save Table A15 as .pdf
pdf("tableA15.pdf", height=14, width=8.5)
grid.table(full_results_table)
dev.off()


